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Several brain diseases are characterized by abnormal neuronal synchronization. 
Desynchronization of abnormal neural synchrony is theoretically compelling because of 
the complex dynamical mechanisms involved. We here present a novel type of coordinated 
reset (CR) stimulation. CR means to deliver phase resetting stimuli at different neuronal 
sub-populations sequentially, i.e., at times equidistantly distributed in a stimulation cycle. 
This uniform timing pattern seems to be intuitive and actually applies to the neural 
network models used for the study of CR so far. CR resets the population to an 
unstable cluster state from where it passes through a desynchronized transient, eventually 
resynchronizing if left unperturbed. In contrast, we show that the optimal stimulation times 
are non-uniform. Using the model of weakly pulse-coupled neurons with phase response 
curves, we provide an approach that enables to determine optimal stimulation timing 
patterns that substantially maximize the desynchronized transient time following the 
application of CR stimulation. This approach includes an optimization search for clusters 
in a low-dimensional pulse coupled map. As a consequence, model-specific non-uniformly 
spaced cluster states cause considerably longer desynchronization transients. Intriguingly, 
such a desynchronization boost with non-uniform CR stimulation can already be achieved 
by only slight modifications of the uniform CR timing pattern. Our results suggest that the 
non-uniformness of the stimulation times can be a medically valuable parameter in the 
calibration procedure for CR stimulation, where the latter has successfully been used in 
clinical and pre-clinical studies for the treatment of Parkinson's disease and tinnitus. 

Keywords: desynchronization, pulse coupled neurons, coordinated reset stimulation, phase response curve, 
stimulation timing, cluster states 



1. INTRODUCTION 

Pathological neuronal synchronization is a hallmark of several 
neurological disorders like Parkinson's disease (PD), essential 
tremor, tinnitus, or epilepsy (Lenz et al., 1994; Nini et al., 1995; 
Mormann et al, 2000; Weisz et al, 2005, 2007; Kane et al., 2009; 
Schnitzler et al., 2009; Roberts et al, 2010), whereas the neu- 
ronal firing is uncorrelated in the normal state (Nini et al, 1995; 
Wilson et al., 2004) such that the abnormal synchronization is 
associated with pathology and symptoms (Levy et al, 2000). 
The standard therapy for medically refractory PD patients is 
electrical deep brain stimulation (DBS), where a high-frequency 
(HF > 100 Hz) electrical pulse train is administered to target 
brain areas via chronically implanted depth electrodes (Benabid 
et al., 1991). HF DBS is found to significantly alter the neu- 
ronal activity of the stimulated neurons: The neuronal firing 
can be suppressed in the vicinity of the stimulation electrode, 
whereas the neurons are overactivated in the output structures 
of the stimulated neuronal population such that the patho- 
logical firing pattern is changed, see (Deniau et al., 2010) for 
review. 

The mechanism of HF DBS is not fully understood. A 
modeling study by Wilson et al. (2011) suggests that the HF 
periodic DBS may induce a chaotic desynchronization, while 



a desynchronizing impact of a periodic forcing on synchro- 
nized populations seems to be a rather general phenomenon 
(Popovych and Tass, 2011). Also, as shown computationally, 
the effect of HF DBS strongly depends on the target struc- 
tures stimulated (Hauptmann and Tass, 2007): Delivering HF 
DBS (nearly) exclusively to excitatory target structures may 
cause a desynchronization, whereas a stronger involvement of 
inhibitory target structures typically causes a pronounced inhi- 
bition. In some patients, however, HF DBS may be ineffective 
or cause side effects (Limousin et al., 1999; Kumar et al., 2003; 
Volkmann, 2004; Freund, 2005; Rodriguez-Oroz et al, 2005). 
Accordingly, along the lines of a model-based approach (Tass, 
1999) novel stimulation techniques have been developed (Tass, 
2001a, 2003a,b; Rosenblum and Pikovsky, 2004; Hauptmann 
et al, 2005; Popovych et al, 2005, 2006; Pyragas et al, 2007; 
Popovych and Tass, 2010), which selectively counteract the patho- 
logical synchronization and restore uncorrelated neuronal firing. 
They are based on either phase resetting or feedback approaches, 
where from a theoretical standpoint the latter might have a great 
potential in controlling pathological synchronization, but so far 
have not been technically realized. 

Other methods suggest to stimulate a single oscillator in the 
network (Nabi and Moehlis, 2011), drive the neurons into a 
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phaseless set in order to achieve desynchronization (Danzl et al., 
2009), or focused on the optimization of the standard HF DBS via 
a closed-loop stimulation setup (Feng et al., 2007a,b). In monkeys 
rendered parkinsonian with the neurotoxin MPTP Rosin et al. 
(2011) studied closed-loop DBS under acute conditions. To this 
end, they delivered a short train (comprising 7 pulses at 130 Hz) 
through a pair of electrodes located in the Globus pallidus inter- 
nus (GPi) at a predetermined, fixed latency (80 ms) following 
each action potential recorded through an electrode placed in the 
primary motor cortex (Ml). This type of stimulation caused a 
strong decrease of the firing rate of the pallidal neurons together 
with a pronounced decrease of the oscillatory neuronal activity 
at tremor frequency (4-7 Hz) and at double tremor frequency 
(9-15 Hz) along with an amelioration of the MPTP-induced aki- 
nesia. After cessation of this type of closed-loop DBS the initial 
firing pattern reverted back, i.e., pallidal firing rate and palli- 
dal oscillatory activity attained pre-stimulus levels (Rosin et al., 
2011). In contrast, standard continuous 130 Hz DBS caused a 
less pronounced decrease of the pallidal firing rate, the oscilla- 
tory neuronal activity and the amelioration of the akinesia (Rosin 
etal, 2011). 

In contrast to the closed-loop DBS tested by Rosin et al. (20 11), 
CR stimulation can be performed in a closed-loop or an open- 
loop mode (Tass, 2003a,b). For instance, an adaptation of the 
stimulation frequency to the dominant frequency of the patho- 
logical neuronal synchronized collective oscillation can increase 
its efficacy (Tass, 2003a,b), see also (Tass et al., 2009). However, 
CR stimulation is robust with respect to variations of both stimu- 
lation and model parameters as follows from both computational 
as well as pre-clinical and clinical studies (Tass, 2003b; Tass et al., 
2012a,b). More importantly, the goal of the CR approach is fun- 
damentally different. CR stimulation does not aim at a decrease of 
firing rates and/or an abolishment of oscillatory neuronal activity. 
Rather, CR stimulation aims at specifically counteracting patho- 
logical synchrony by desynchronization (Tass, 2003a,b). This is 
because neurons have to be active in order to unlearn their patho- 
logical synaptic connectivity. In this way sustained long-lasting 
desynchronization is induced, and, as predicted computation- 
ally (Tass and Majtanik, 2006; Hauptmann and Tass, 2007; Tass 
and Popovych, 2012), therapeutic effects were observed in rat 
hippocampal slice experiments in the context of epilepsy (Tass 
et al., 2009) as well as in a clinical proof of concept study in tin- 
nitus patients treated with acoustic CR stimulation (Tass et al., 
2012a), where both pathological neuronal synchrony (Adamchic 
et al., 2013) and pathological effective connectivity (Silchenko 
et al., 2013) considerably decreased. In addition, in parkinsonian 
(MPTP) monkeys it was shown that unilateral CR stimulation 
delivered to the subthalamic nucleus (STN) for only 2 h per day 
during 5 days leads to significant and sustained therapeutic after- 
effects for at least 30 days, while standard 130 Hz DBS has no 
aftereffects (Tass et al., 2012b). Another motivation why CR stim- 
ulation is delivered at frequencies similar to the pathological oscil- 
latory frequency, is that in this case the desynchronizing effect is 
achieved at favorably small stimulation intensities (Tass, 2003a,b). 
In fact, as shown computationally, CR stimulation is able to 
strongly alter neuronal firing rates if delivered at frequencies sub- 
stantially different to the dominant frequency of the stimulated 



neuronal population (Lysyansky et al., 2011b). For instance, CR 
stimulation may effectively activate hypo- or inactive neuronal 
populations without inducing neuronal synchrony (Lysyansky 
etal, 2011b). 

The CR stimulation (Tass, 2003a,b) is based on the phase 
reset of oscillatory neuronal activity and has a broad applicabil- 
ity, since the phase reset is a universal phenomenon and can be 
achieved for a variety of stimulation setups and conditions, see, 
e.g., Refs. (Winfree, 1977; Paydarfar and Eldridge, 1987; Brandt, 
1997; Makeiget al, 2002; Neiman et al, 2007; Thorne et al, 2011). 
According to its stimulation protocol, CR stimulation counteracts 
synchronization in the neuronal target population by dividing 
the entire population into several sub-populations where the 
phases of the neuronal oscillators within each sub-population are 
reset by the stimulation sequentially, i.e., in a timely coordinated 
manner. In this way, the collective neuronal oscillations in the 
sub-populations get phase shifted with respect to each other, and 
the total synchronization is replaced by, e.g., a cluster state (Tass, 
2003a,b; Lysyansky et al., 201 la). Due to the pathologically strong 
synaptic connectivity, the entire target population runs from the 
cluster state through a transient characterized by pronounced 
desynchronization and finally resynchronizes if left unperturbed. 
Accordingly, to keep the neuronal ensemble in a desynchronized 
state, CR stimuli are delivered intermittently (Tass, 2003a,b), for 
instance, by applying CR in an m:n ON-OFF mode, where a 
few, say m, cycles with CR are optimally followed by a few, say 
«, cycles without any stimulation (Lysyansky et al., 2011a), for 
example, with m = 3 and n = 2 (Tass et al., 2012a,b). Such a stim- 
ulation protocol has computationally been found to be effective 
in inducing transient desynchronization in stimulated neuronal 
ensembles (Tass, 2003a,b; Lysyansky et al., 201 la). 

The post-stimulus transient, where the stimulation-free neu- 
rons undergo an unperturbed desynchronized dynamics, plays an 
important role for the emergence of long-lasting effects of CR 
stimulation. In computational models taking into account the 
adaptive synapses governed by spike timing-dependent plastic- 
ity (STDP) (Gerstner et al, 1996; Markram et al, 1997; Feldman, 
2000; Wittenberg and Wang, 2006; Caporale and Dan, 2008), it 
has been shown that CR stimulation can lead to a reduction of 
the mean synaptic weight and, in turn, shift the network to a state 
characterized by desynchronized activity and weak connectivity 
(Tass and Majtanik, 2006; Hauptmann and Tass, 2007) which per- 
sists after the stimulation is switched off. Modeling shows that CR 
stimulation is effective for a number of stimulation setups, in par- 
ticular, for direct somatic stimulation as well as for excitatory or 
inhibitory synaptically-meditated stimulation which corresponds 
to stimulation of afferent or efferent fibers (Popovych and Tass, 
2012). This is particularly important since it has been shown that 
stimulation of fibers projecting to the STN appears to be respon- 
sible for the therapeutic effect of HF DBS delivered through STN 
electrodes (Gradinaru et al, 2009). 

In order to further optimize the therapeutic benefit of CR 
stimulation, in this paper we investigate the impact of the stimula- 
tion parameters and the stimulation protocol on the stimulation- 
induced desynchronization. In particular, we focus on how the 
timing of the phase resets of the neuronal sub-populations influ- 
ences the quality of the stimulation-induced cluster state and the 
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post-stimulation transient. We found that appropriately adapted 
non-uniform stimulation onsets for the different stimulation sites 
can divide the phases of the stimulated neurons in such a way that 
the desynchronized post-stimulation transient gets significantly 
prolonged, until the population eventually resynchronizes again. 
To confine the complexity of our analysis, we study phase mod- 
els without STDP. Phase resetting can be incorporated in such 
models in a natural way. We here set out to determine an opti- 
mal pattern of phase resets for CR stimulation. Put otherwise, we 
address the question of how to optimally choose the stimulation 
onsets for the single stimulation sites in CR. 

In weakly coupled networks of oscillators, the technique of 
averaging is often applied to obtain a coupling which involves 
only the relative phase differences of the interacting oscillators 
(Ermentrout and Kopell, 1991; Swift et al., 1992; Hoppensteadt 
and Izhikevich, 1997). In globally coupled systems of identical 
neurons, these averaged phase models always possess symmetric 
cluster states, i.e., states, in which m clusters of equal size exist 
with a phase distance of 2tc /m between neighbors (Okuda, 1993). 
Hence, a natural answer to the question above is to choose the 
stimulation times such that the phases get uniformly distributed 
over one period, independently on the type of neurons. In this 
case target patterns of the stimulation are the symmetric clus- 
ter states. However, in non-averaged models, the importance of 
the coupling between the single neurons becomes apparent. It 
plays an important role in determining the exact way of how 
the stimulation should be applied to cause a longer transient. In 
this work, we use systems of globally pulse-coupled phase oscil- 
lators for modeling the dynamics of a neuronal population. In 
particular, the symmetric cluster states disappear generically, and 
non-symmetric cluster states become possible candidates as stim- 
ulation target states. We propose a method for computing the 
stimulation times, which resets the system to a suitable cluster 
state. The timing points of the applied stimuli in these cases 
are non-uniformly spaced. The desynchronizing post-stimulation 
transient after such a stimulation turns out to be longer than 
the corresponding post-stimulation transient after a uniform CR 
stimulation of the same system. 

The paper is organized as follows. In section 2.1 we intro- 
duce the globally pulse-coupled model that is used to describe 
the collective dynamics of the neurons. In section 2.2 we study 
the relevant dynamical properties of the model, i.e., the appear- 
ance of synchronization, symmetric clusters and splay states. 
CR stimulation is introduced to the model in section 2.3, and 
we derive how one should apply CR stimuli to obtain longer 
post-stimulation desynchronization transients in section 2.4. 
The theoretical analysis is illustrated by numerical simulations 
in section 3. In particular, the robustness of the results to 
the variation of the stimulation parameters (stimulation inten- 
sity and electrode activation time) is studied in section 3.2. 
The effects of inhomogeneous frequencies is studied in sec- 
tion 3.3. 

2. MATERIALS AND METHODS 

2.1. PULSE-COUPLED PHASE OSCILLATORS 

Phase models play a key role in describing the individual 
dynamics of single oscillators, e.g., oscillatory neurons, see e.g., 



(Kuramoto, 1984). In particular, a stable periodic dynamics can 
be modeled by a simple equation for the periodic motion of the 
phase with frequency co: tp(f) = cp (0) + oof, where (p is consid- 
ered modulo 2tt. The direction of the phase is neutrally stable. 
Therefore, a sufficiently weak temporary perturbation, which 
does not move the original system far away from the correspond- 
ing limit cycle, persists in the phase for all times, while all its other 
effects die out exponentially fast due to the stability of the limit 
cycle which corresponds to the periodic motion. In coupled sys- 
tems, weak interactions can be conceived as perturbations, and 
the phase reduction can be applied as well (Kuramoto, 1984; 
Hansel et al., 1993; Hoppensteadt and Izhikevich, 1997; Brown 
et al., 2004). In fact, phase models are particularly important 
for studying network dynamics, because many types of synchro- 
nization, which are of interest in such models, depend on the 
relative phases of the units combining the network (Pikovsky 
et al, 2001). 

The effect of perturbations is incorporated into the phase 
equations by the phase response curve (PRC) (Guckenheimer, 
1975; Kuramoto, 1984; Ermentrout, 1996; Winfree, 2001) (see 
Figure 1). It measures the response of the individual neu- 
ron to weak stimuli. We consider the case when it admits 
a representation in terms of a smooth scalar function Z(cp) 
of the phase. For example, applying a weak current 1(f) to 
a neuron with PRC Z((p) changes its phase dynamics from 
ip(f) = co to 

<p(f) = oo + I(f)-Z(cp(f)), 

which describes the phase dynamics in the weak stimulation limit 
(Kuramoto, 1984). If the perturbation is pulse-like, e.g., a brief 
electrical stimulation or a synaptic input from another neuron 
(if synapses are fast), it may be approximated as an instanta- 
neous input which resets the neuron's phase at time t = frj of the 
incoming pulse as following 

<p » <p ( f o + ) = <p for) + 1 ■ z Op (<o~)) > 

where <p(f^) = linittfo <p(f) and cp (fj~ ) = lim ( j, to cp(f). We can 
formally write the pulse-like perturbation as 7(f) = I ■ 8 (f — fn), 
using the Dirac delta-function (Goel and Ermentrout, 2002). In 
what follows, we study a system of N identical phase oscillators 




0 (p 2it 0 <p 271 



FIGURE 1 | Phase response curves. (A) PRC Z(cp) = Z H (<p) = - sin(ip) of 
an oscillator close to a bifurcation of Hopf- or Bautin-type (Brown et al., 
2004). (B) PRC Z(cp) = Z M L(ip) of a particular Morris-Lecar type model 
(Ermentrout, 1996; Sato et al., 2011). 
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which are globally pulse-coupled with weight I = (Goel and 
Ermentrout, 2002) 



N 



(i) 



k=l I 



where tk,t are the times where the A:-th neuron spikes. For conve- 
nience we assume that a single neuron emits a spike if its phase 
crosses zero (mod lit). In such a way, if the phase % of neu- 
ron k passes through the spike, cpjt (t^^j = 2jt, all neurons q>j are 

reset to (p,- (^li) = ^ yPj ( f i<^)) wnere me resetting function (x 
is defined as 



M-(<P) 



<P + ^ -Z(cp). 



(2) 



In the case that n neurons simultaneously spike at time t = f*, 
the reset value of neuron (p<- is taken to be cp, = u," (cp,- (i*r)), 
where (Ji"(<p) = |x (|x (. . . u,(cp) . . .)) denotes the n-fold super- 
position of the resetting function. This mapping is preferable 
to the choice of ip ; - = cp ; - (t~) + '-f^Z (tpj (**")) because it 
assures a continuous dependence on the initial conditions. More 
specifically, the effect of the spiking of a cluster of n neu- 
rons changes continuously as the neurons split from the cluster, 
i.e., u, (ti + \l (X2 + ]x (. . . jJt(<P) • • ■))) M-"(<P) as the interspike 
intervals T; — > 0. Note that the resetting function for the clus- 
ter spikes can also be measured or computed directly in order to 
achieve more realistic modeling (Achuthan and Canavier, 2009). 
Applied to our case it would mean using a given measured func- 
tion u,(cp, n) instead of u,"(cp) in the case when an M-cluster spikes. 
However, we restrict our analysis to the function u"(<P) defined as 
a superposition of individual resetting functions. 

Note, that for sufficiently small values of x/N the resetting 
function \i is monotone. This ensures the preservation of the 
ordering of the phases. We assume that N is sufficiently large 
and this property holds. For some systems (Goel and Ermentrout, 
2002; Brown et al., 2004; Stiefel et al, 2009), it turned out, that the 
PRC has a small value at the spike moment Z(0) » Z(2it) ~ 0. 
For simplicity, we assume that Z(0) = Z(2tc) = 0 holds for our 
model. Let us shortly explain why this is a reasonable approxima- 
tion for weakly coupled spiking systems. Firstly, the modulus of 
the PRC must be roughly proportional to the density of isochrons 
(Guckenheimer, 1975; Winfree, 2001) at the corresponding point 
of the limit cycle of the full system. This density on the other hand 
is inversely proportional to the modulus of the vector field, which 
is large at the spiking point. Therefore, the modulus of the PRC 
has to be small at the spiking point. 

Note, even though this work is focused on systems with pulse 
coupling (1), the main qualitative message about non-uniform 
CR stimulation and non-uniform positions of clusters still holds 
for systems of the form 



(3) 



with a smooth, periodic coupling function G, which was 
proposed in (Winfree, 2001). Systems (1) and (3) were less 
extensively studied than their averaged versions, which take the 
form 



(4) 



where H(<p) = (2jt) _1 / 0 "2(i|f)G(«p + i|f),Ji|;, see, for exam- 
ple, (Ermentrout and Kopell, 1991; van Vreeswijk et al., 1994; 
Daido, 1997; Kuramoto, 1997). As a result of the averaging, 
the stability properties of the corresponding solutions of (1) 
[or (3)] and (4) may differ at the order of O (x 2 ). This is the 
same magnitude as of the errors made by the phase reduc- 
tion, and, thus, studying the averaged system suggests itself as 
a simpler and, presumably, equivalent task. In the next sec- 
tion we show that an important genericity of stationary solu- 
tions with distributed phases is overlooked by this choice. In 
fact, the homogeneous stationary solutions or symmetric clus- 
ters in the averaged system (4) correspond to some other, 
generically non-homogeneous solutions of the original system, 
whose shape may differ at the order of O (h), see, for instance, 
Equations (12) and (13) below. A precise targeting of these 
solutions by CR can essentially contribute to the efficacy of 
the desynchronization technique. Since systems (1) and (3) 
admit non-homogeneous stationary solutions with distributed 
phases, it is of particular interest to study them in the context 
ofCR 

2.2. SYNCHRONIZATION, CLUSTERS. AND STATIONARY SOLUTIONS 

In this section we review the dynamical properties of the 
stimulation-free population (1), which are relevant in the con- 
text of desynchronization. In particular, we study the appear- 
ance and stability of a synchronized state, cluster states, as 
well as splay states. We pay special attention to the differences 
between the mentioned dynamical states of models (1), (3), 
and (4). 

2.2. 1. Stability of synchronized solutions 

In each of the systems (1), (3) and (4), there exists an in-phase 
synchronous solution where all neurons are perfectly synchro- 
nized 



cpi(0 



<Pn(0- 



The conditions for the stability of the synchronous state are 
well known for all of the above systems (Goel and Ermentrout, 
2002). In particular, for the pulse-coupled system (1), in-phase 
synchronization is linearly stable iff 



xZ'(0) < 0. 
For (3) the corresponding condition reads 



Jo 



G((p)Z'((p) 
oj + xZ((p)G((p) 



ci(p < 0. 



(5) 



(6) 
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For the averaged system (4), the stability condition of the syn- 
chronized state is 

- xH'(O) = — / G(<p)Z'(<p)d<p < 0. (7) 
2tt Jo 

A comparison of conditions (5)-(7) leads to the following rela- 
tionships: 

- Vanishing coupling: Conditions (7) and (6) for the averaged 
and non-averaged systems differ only in the second order of 
x, respectively. Therefore, they coincide in the limit of small 
coupling x — > 0, if H' (0)^0. 

- Smooth, pulsatile coupling: If the coupling function G(cp) is 
pulse-like, i.e., it is positive and concentrated at tp = 0, then 



and for the averaged system (4) one gets 



H((p) % GZ(-cp), G 



2jt Jo 



and the condition for synchronization of the averaged system 
(7) coincides with the condition for the non-averaged pulse- 
coupled system (5) provided Z'(0) ^ 0. 

Although in this work we will focus on smooth PRCs, it is 
worth to mention another synchronization effect, which may 
generically occur for PRCs, which have discontinuous deriva- 
tives as a phase of zero is approached from the left and right 
(Lucken and Yanchuk, 2012): though the synchronous state is 
locally unstable in this case, the first order parameter R\(t) = 
I T! 12k ex P ( zc PJc(0)| may still approach its maximal value R\ = 
1 due to the appearance of structurally and dynamically stable 
homoclinic connections to the synchronous state. 

2.2.2. Splay states and stationary solutions 

Splay states are periodic solutions of system ( 1 ) , in which all oscil- 
lators are spread in a way that the time differences between the 
subsequent spikes tk+i t i — tk,i are always the same (Swift et al., 
1992; Zillmer et al., 2007). Note that the term "splay state" can 
also be used differently and may, more generally, refer to any 
state featuring phases which are spread over the periodic interval 
[0, 2jt] (Achuthan and Canavier, 2009). Splay states are impor- 
tant for desynchronization issues, since they possess small order 
parameters R n =\jjJ2k ex P ( in Vk(t))\. 

To study splay states in large systems it is useful to consider 
an equation for the phase distribution density p(f, cp), since its 
stationary solution approximates the distribution of the phases of 
splay states in the limit of large N. For the pulse-coupled system 
(1) the dynamics of the phase distribution density is governed by 
the following continuity equation (Ernst et al., 1995; Brown et al., 
2004) 

3 f p(f, cp) = -3 V (p (t, cp) (oo + xZ(cp)p (t, 0))) . (8) 
Its equivalent for the smoothly-coupled model (3) is 



d t p (t, cp) 



(f) (t, <P) (u> + xZ(cp) 



G (\|0 p (t, \|/) d^ 



fin 



oj + x J H(i\r - cp) p(t,\\i)dty 

(10) 



Solving (8) for stationary solutions p(t, cp) = p 5 (cp) and taking 
into account the normalization f Q 71 p s (ty)dty = 1, we obtain 



2jt 



oj + xZ(cp)p s (0) 
oo + xZ 010 p s (0) 



(11) 



Here p s (0) is uniquely determined by the implicit equation 
obtained from (11) by inserting cp = 0 (see Appendix A.l for 
details). Thus, (11) describes a unique stationary solution of (8). 
For small x, this solution can be approximated as 



Ps(<P) 



Z(<p) 



1 Z 
27t (2jt) / oo 



+ O (x 2 ) , 



(12) 



where Z : 



2n JO 



Z (\|x) ci\|/ is the mean value of the PRC. For 



smoothly-coupled systems (3), one analogously finds that a 
unique stationary solution of (9) exists, if x is not too large (see 
Appendix A.l ). Its first-order expansion in x reads 



Ps(<P) 



1 +x G(Z-Z0p)) +o(x2) 



2jt 



2ttoo 



(13) 



For the averaged model (4), we find that for any value of x the 
constant distribution density p(cp) = is a stationary solution 
of(10). 

As follows from Equations (12) and (13), phase distributions 
of the splay states of the pulse-coupled system (1) as well as the 
non-averaged system (3) deviate from a uniform distribution. For 
small x, the deviations can be estimated in the first order in x by 
(12) and (13), respectively. This is in contrast to splay states of 
the averaged system (4), which are always uniformly distributed. 
Figures 2A,B illustrate non-uniform stationary phase distribu- 
tions in pulse-coupled systems with PRCs Zn(<p) (Figure 1A) 




(9) 



FIGURE 2 | Stationary phase distribution densities. Black solid curves 
show the stationary phase distribution densities p s (ip) (with the scale on 
the left vertical axes) for the pulse-coupled systems (1) with PRCs (A) 
Z(cp) = Z w (<p) (Figure 1A) and (B) Z(cp) = Z ML (ip) (Figure IB) The 

corresponding PRCs are depicted by gray dashed curves and rescaled by 
some constant ratio (with the scale on right vertical axes). Coupling 
strength x = 1 .0. 
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and Zml(<P) (Figure IB), respectively. Figure 3 shows that the 
theoretically obtained stationary phase distribution density p 5 (ip) 
(black curve) convincingly approximates the numerically calcu- 
lated phase distribution histogram (gray bars) of the splay state in 
the pulse-coupled system (1) for large number N of oscillators. 

The stability of splay states as well as cluster states can be stud- 
ied by the set of multipliers X of the corresponding return map, 
which determine the rate, with which small perturbations from 
the considered state are growing with time. In particular, if the 
multiplier with the maximal absolute value is X max > then a generic 
perturbation will grow as \ max with j is the number of spikes 
of a neuron from the ensemble. Therefore, a stable state corre- 
sponds to the case |X max | < 1 and an unstable to |X max | > 1. 
We calculated numerically the maximal multipliers X max of the 
splay state of system (1) for two cases: (i) Z(cp) = Z/j(<p) and (ii) 
Z(cp) = Zml(H>) (Figure 4, blue pluses). We observe that the splay 
states are unstable for all N except N = 3 in case (ii), and their 
multipliers converge to an asymptotic value Xqq > 1 as the num- 
ber of neurons N increases. Note that the splay state may be stable 
for synaptic coupling and type-I PRC, as shown in (Achuthan 
and Canavier, 2009) for a system of a few coupled oscillators. For 
a more detailed analysis of the stability of splay states in large 



pulse-coupled systems, we refer to (Abbott and van Vreeswijk, 
1993; Zillmer et al, 2007; Calamai et al, 2009). 

2.2.3. Symmetric clusters 

In this paper we consider strongly synchronized neuronal ensem- 
bles, where the splay state is expected to be unstable. When an 
external desynchronization technique would be able to move the 
system in a vicinity of such a state, the achieved order param- 
eter would be very low for a relatively long time. In the case, 
when the number of stimulation sites is naturally limited to 
a low number, a natural substitute for the target state of CR 
stimulation is a cluster state. For the sake of simplicity we will 
consider symmetric cluster states consisting of m clusters, each of 
them containing N/m neurons. Within each cluster, the neurons 
are synchronized and have the same phase, whereas the phases 
of different clusters are shifted with respect to each other. For 
systems (1) and (3) there exists a unique stationary, symmet- 
ric m-cluster state, at least for moderate coupling strength. The 
phases \|/ ; - of individual clusters are not equidistantly distributed, 
i.e., in general (f) — \py(t) | ^ 2n/m. In contrast, in the aver- 
aged systems of the form (4), equidistantly distributed m-clusters 
with |i|f; +1 (f) — = 2n/m always exist. In Appendix B, we 

explain how one can determine multipliers of symmetric clus- 
ter states for system (1). We have computed them for the cases 
(i) Z((p) = Z H (y) and (ii) Z(<p) = Z M l(<P) (Figures 4A,B). It is 
convenient to distinguish between tangential and transverse mul- 
tipliers, which correspond to perturbations within the invariant 
cluster space and perturbations which disperse the single clusters 
by destroying the perfect synchrony within the clusters, respec- 
tively. For case (i), and 1 < m < 25, all m-cluster states are 
unstable. In case (ii), there is a single value m = 3 for which the 
cluster state is stable. The tangential multipliers seem to asymp- 
tote towards a limit with increasing m (Figure 4, red circles). If 
the transverse multipliers are smaller than the tangential, as it 
is the case for both PRCs (Figure 4, red markers), one expects 
that the perturbed dynamics stays near the invariant cluster space 
as long as the linear prediction is valid. Detailed analysis of clus- 
ter states for different PRCs has been performed in (Ashwin and 
Swift, 1992; Hansel et al, 1993; Okuda, 1993; Chandrasekaran 
et al, 2011; Lucken and Yanchuk, 2012). 
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FIGURE 3 | Comparison of splay state and stationary phase 
distribution density. Black solid curve depicts the theoretical solution of 
the stationary phase distribution density p s (ip) of (8) for the pulse-coupled 
system (1) with a PRC of Morris-Lecar type Z(<p) = Z ML (<p) (Figure 1B). 
Gray bars illustrate the numerically computed and normalized phase 
distribution histogram of the corresponding splay state for N = 1000 
oscillators. Coupling strength u = 3. 
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FIGURE 4 | Characteristic multipliers of splay states and cluster 
states. Maximal moduli \ m ax of characteristic multipliers of cluster 
states (red "o": tangential, red "A": transverse) and splay states 
(blue "+") versus the number of clusters m, respectively the number 



of neurons N, for (A) PRC Z = Z H and (B) PRC Z=Z ML . Multipliers for 
cluster states have been calculated asymptotically for large N using an 
asymptotic technique described in Appendix B. Coupling strength 
k = 0.5. 
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2.3. MODELING DEEP BRAIN STIMULATION 

Strong enough electrical stimuli or synaptic input can reset the 
phase of a neuron in such a way that its oscillation restarts after 
the stimulation from a definite phase (Winfree, 1977; Best, 1979; 
Tass, 1999; Popovych and Tass, 2012). The general mechanism 
of phase resetting can be understood as follows. In a stimulated 
neuron, a stable steady (e.g., hyperpolarized) state appears, which 
is approached during the stimulation. When the stimulation ter- 
minates, the steady state disappears, and the system relaxes back 
to the limit cycle with an asymptotic phase (p s determined by 
the isochron on which the stimulation specific steady state was 
located. However, for the phase resetting procedure which we pro- 
pose, the exact value of cp s is not essential. Therefore we include 
the stimulus in the simplest way which provides the model with 
a qualitatively adequate stimulus response. The phase model with 
the incorporated stimulus reads 



\ k=l I I 



(14) 



where lj(t) corresponds to the stimulus applied to the;'-th neuron. 
We assume, that there are m stimulation sites, each stimulating a 
distinct group of N/m neurons (Figure 5). These stimulation sites 
can either be active or not, that means 



Ij(f) 



I, fort eSj = \Jilj,i 



0, else, 



(15) 



where X;j are stimulation time intervals and I is the stimulation 
intensity. For two neurons (p; and <p* from the same group we have 
Sj = S k and I } (t)=I k (t). 

The dynamics of uncoupled neurons under stimulation is 
described as 



< + lZ(ipj). 



(16) 



If IZ(<p) < 0 for some (p, and if the intensity I is of sufficiently 
large magnitude, there appears a pair of fixed points, stable <p s 
and unstable cp„, satisfying oj + IZ(cp) = 0. If only one such pair 
exists, the neuron will approach the stable fixed point <p s after 
some time of stimulation and stay there until the stimulation ter- 
minates (see Figure 6). In such a situation, we call (p s the resetting 
point. 

The stimulation described by (14) aims at establishing a dis- 
tribution of phases of the neuronal ensemble that prolongs the 




FIGURE 5 | Schematic illustration of the CR stimulation setup, m = 4 

stimulation sites (arrows with different filling patterns) affect the 
corresponding distinct sub-populations of N/m neurons (circles with the 
same filling pattern). 



post-stimulus transient as much as possible before the ensem- 
ble synchronizes again. In principle, the strategy is to establish 
a state as close as possible to some stationary desynchronized 
state. With a given number m of stimulation sites which influ- 
ence equally large groups of N/m neurons, the target states for 
the control are restricted to m-cluster configurations with clus- 
ters of equal size. We call a target pattern the state, which is 
intended to be realized at the end of the stimulation. A series of 
successive activations and deactivations of the stimulation sites 
is called stimulation sequence, and a time interval during which 
the resetting stimuli are delivered at all tn sites is called stimula- 
tion cycle. For the averaged system (4), the target pattern consists 
of equidistant clusters = + (p s , 1 < j < m, such that the 
last stimulation-induced cluster is located at the resetting point 
\[r m = (p s at the end of the stimulation cycle. For systems of type 
(1) or (3), the target pattern is a cluster state = 7?/+ <Ps + 
0(x), which is in general not equidistantly distributed. Among 
all possible stimulation sequences, we restrict our considerations 
to those where each stimulation site is activated once per stimu- 
lation cycle. The activity of the j-th stimulation site is confined 
to the time interval Sj = [fj, tj + t], where t\, . . . , t m are the 
onset times within the stimulation cycle, and x is the stimulation 
duration, which is the same for all stimulation sites. In practice, 
stimulation sequences have to be administered repeatedly after 
the system recovers to some undesired level of synchronization 
(Tass, 2003a). 

2.4. STIMULATION-INDUCED STATIONARY M-CLUSTER STATES 

Since the stimulation target pattern has to be established as pre- 
cisely as possible, one has to take into account the influences of 
the coupling among neurons on the stimulation-induced pat- 
tern. We now describe how appropriate stimulation sequences 
can be found when we restrict ourselves to the stimulation timing 
Sj = [tj, tj + x], j = 1, . . . , m, mentioned above. The long-term 
desynchronizing effects of such a stimulation sequence will be dis- 
cussed in the next section 3. For the brevity of notations, let us 
introduce the resetting map of the stimulated system (14) by 

$ : R'" x [0, 2ti] n -> [0,2jt] n , 
(t; cp) <D(t; <p), 




FIGURE 6 I Phase dynamics of the uncoupled neurons (16) during 
stimulation. (p s and <p u denote stable and unstable fixed points, 
respectively. The arrows indicate the direction of convergence of the phase 
to the stable fixed point <p s . PRCs and stimulation intensities (A) 
Z(<p) = Z w (cp), /= 10 and (B) Z(<p) = ZmlM, / = -10. 
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where (p e [0, 2n] N is the state of the system at the onset of the 

def 

stimulation at t = t on = min{ti, . . . , f m }, and t = (ti, . . . , t m ) is 
the vector of the onset times of the stimulation sites. 0(t; (p) 
describes the JV-dimensional state of the system at the end of 
the stimulation. If the duration x and magnitude of the stimu- 
lation intensity I are large enough (see section 3.2), each neu- 
ron cpi- of the j-th group is reset to the collective cluster phase 
cpjt (tj + r) ~ (p s at the offset time fc + x of the j-th stimulation 
site and continues to evolve in a cluster of a common phase 
which we denote by \[fj (f; cp). The dependence of on initial 
conditions (p follows from the fact that this cluster may still be 
influenced by the other clusters and neurons. In practice, the 
state (p before stimulation is unknown. To ensure full control 
over the resulting cluster state *(t; cp), it should be indepen- 
dent of the initial state (p. To this aim we have to assume that 
the resetting mechanism results in an accurate reset of the stim- 
ulated neurons to a determined phase (p = <p s . Furthermore, we 
have to ensure that the reset neurons are not affected by neurons, 
which have not yet been reset in the current stimulation cycle. 
Both conditions can be fulfilled by choosing a sufficiently large 
duration x and intensity I for the reset. Thus, we can approx- 
imately identify the state <J>(t; cp) with the lower dimensional 
cluster state * (t) = (f off ) , . . . , ty m (t 0 s)), where f off = x + 
max{fi, . . . ,t m ] denotes the offset time of the entire stimula- 
tion sequence, and moreover, * (t) is independent on (p. Now, 
let = (\|/*, . . . , 4>f„) denote the phases of the stationary clus- 
ter state which serves as the stimulation target pattern and can be 
obtained as described in Appendix B. Then, the problem is to find 
a solution t of 



*(t) = **. 



(17) 



In this study we do not aim to provide a general algorithm to solve 
this equation. Moreover, due to discontinuities that are caused by 
the pulse-like interactions of the ensemble, the function * (t) is 
only piecewise smooth (see Figure 8). We approach the solution 
of (17) numerically by starting from the uniformly distributed 
stimulation sequence to as an initial guess. Then we apply the 
minimization Nelder-Mead simplex search algorithm (Lagarias 
et al., 1998), which is implemented in the MATLAB function 



Table 1 | Example of target pattern (stationary cluster state) 
** — (i|r*, . . . , i|r*) and stimulation sequence t = (tj, .. . , tj). 
Parameters: m = 4, Z = Zml, w = 1, x = 0.5, x = 2n, I = — 10, and 
ton = 0. 



1.169690 
2.703008 
4.284900 
5.918402 



0 

4.749867 
3.215394 
1.633501 



At,. 



1.534474 
1.581893 
1.633501 



The time elapsed between t* and the preceding stimulation onset is denoted 
by A t*. Note the non-uniformity of the stimulation. See also the corresponding 
illustration in Figure 7 



fminsearch, to minimize [|* (t) — **||. Table 1 and Figure 7 
illustrate an example of the computed stimulation sequence for 
the case of four stimulation sites m = 4 and a Morris-Lecar 
type of PRC Z(<p) = Zml(<P) (see Figure 1). For the parame- 
ter values given in caption to Table 1, the target pattern has 
been computed as well as the stimulation sequence t. The opti- 
mal stimulation sequence deviates by up to ~ 4% from the 
uniformly distributed one where Af* is the same for all = 

2, . . . , m. Figure 7 illustrates the corresponding switching times 
of the stimulation contacts. It also shows the spiking times of 
the obtained clusters after the stimulation. In order to find the 
stimulation sequence, the discontinuous system (17) has been 
solved. Figure 8 illustrates the emerging types of discontinuities 
by plotting ||* (f*, f|, t%, t 4 ) - **| versus f 4 for fixed t*, t*, 
and t%. 

3. RESULTS 

3.1. ADVANTAGES TO UNIFORM STIMULATION 

Equidistant clusters are stationary if the coupling depends only 
on the phase differences as in system (4). In the non-averaged 
systems (1) or (3), the phases of stationary clusters are dis- 
tributed non-uniformly in [0, 2it], and the resetting technique 
described in section 2.4 is expected to yield longer post-stimulus 
transients. The results of numerical simulations of systems (1) 




FIGURE 7 | Illustration of the stimulation sequence, which leads to the 
stationary 4-cluster state. Black horizontal bars indicate time intervals 
when a corresponding stimulation site is active, circles mark the subsequent 
spike times of the established clusters. Parameter values as in Table 1 . 




FIGURE 8 | Illustration of the discontinuity of fc, i-> 4> ((*, tj, (J, f 4 ). 



I*(«f 



, ti) — ** || is plotted versus f 4 for fixed f*, fj, and t^ 



Parameters as in Table 1. The discontinuity occurs at such a value of 
U = td that leads to \|r 4 fc) = 2it, i.e., the onset of the post-stimulation 
spike of the cluster j = 4 just coincides with t e ff, see Figure 7. If £j > t d , 
the impact of the spike of cluster i|/ 4 on the cluster \|/2 has to be taken into 
account when calculating the resulting cluster positions. 
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FIGURE 9 | Advantages of the stimulation-induced target patterns 
consisting of stationary cluster states [plots (B), (D)] compared to those 
induced by conventional CR stimulation with equidistant stimulation 
times [plots (A), (C)]. The time courses of the first order parameter f?i (f) 
(red curves) and the fourth order parameter f?4(f) (blue curves) of the 
neuronal ensemble (14) controlled by CR stimulation are shown for two PRCs 



and stimulation intensities: (A,B) Z(ip) = ZmlM (Figure 1), / = —10 and (CD) 
Z(ip) = Ztf(<p) (Figure 1), / = 10. For each setup, one single CR stimulation 
sequence is administered at f = ton = 250 with duration x = 10 of single site 
activation. Number of the stimulation sites m = 4, system's size N = 240, 
natural frequency w = 1, and the coupling strength x = 0.5. The initial phases 
at f = 0 were randomly drawn from a uniform distribution on [0, 2it]. 



with PRCs Z = Zml and Z = Zh (see Figure 1) are presented 
in Figure 9 where plots (A,C) illustrate the effect of the uni- 
form CR stimulation and (B,D) are related to the non-uniform 
CR stimulation. Time courses of the first order parameter R\ (f) 
(red curves) and the fourth order parameter R 4 (t) (blue curves), 
defined as R n = | ^ ^ k exp (i«(p,t(f))| € [0, 1], are shown. Large 
values of the first order parameter are indicative of an in-phase 
synchronization. On the other hand, for approximately equidis- 
tantly distributed m-cluster states high values of the m— th order 
parameter are combined with low values of all order param- 
eters with lower indices. We use these properties to detect 
synchronization and the discussed slightly non-uniform cluster 
states. 

All simulations in Figure 9 are started at t = 0 with the neu- 
rons' phases randomly distributed in [0, 2n], For both PRCs, 
Z = Zml and Z = Zh, we observe a steady increase of Ri (f) (red 
curves) in the pre-stimulation epoch t < f on , which indicates the 
onset of in-phase synchronization of the entire ensemble. This 
process is significantly faster in the case Z = Zh (Figures 9C,D) 
than for Z = Z M l (Figures 9A,B) due to Z' H (0) « Z ML (0) < 0 
which strengthens the linear stability of the synchronous solu- 
tion in the case Z = Zh (see section 2.2). When CR stimulation 
is turned on at f on = 250, the stimulated phases are successively 



caught at (p s and released when the corresponding stimula- 
tion site is deactivated. At the end of the stimulation cycle at 
f D ff ~ 266.75, four clusters are established, which are well dis- 
tributed in [0, 2tt]. This leads to a low value of J^i (f Q ff) 0 
(A: Ri (t of f) % 0.027, B: Ri (f off ) « 0.020, C: Ri (t oS ) % 0.000, D: 
#i (foff) ~ 0.044) and a high value of R 4 (f off ) « 1 (A: R 4 (f off ) » 
0.987, B: R 4 (t oS ) « 0.994, C: R 4 (f 0 ff) ~ 1-000, D: R 4 (f off ) » 
0.969). All simulations show a post-stimulation transient before 
the system resynchronizes again and the first order parameter 
Ri(t) approaches unity. For both considered PRCs, the advan- 
tages of the method proposed in section 2.4 to establish non- 
equidistant clusters is substantial. For the uniform CR stimulation 
(Figures 9A,C) the post-stimulation desynchronization transient 
is of approximately the same duration as the initial transient in the 
pre-stimulation epoch, when starting from a random distribu- 
tion of the phases. The post-stimulation transient is significantly 
prolonged by the non-uniform CR stimulation in both cases, 
for Z = Zml (Figure 9B) by doubling the transient duration and 
for Z = Zh (Figure 9D) by tripling the duration of the desyn- 
chronization transient. Note that small-scale oscillations of the 
order parameters originate from the discontinuities of the sys- 
tem's trajectory, which occur whenever a cluster crosses the firing 
threshold. 
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3.2. EFFECTS OF DIFFERENT STIMULATION INTENSITIES AND 
DURATION 

An important question is how the stimulation intensity I 
and the stimulation duration x influence the desynchroniza- 
tion transient. Figure 10 shows results of numerical simula- 
tions for the PRC Z = Zml- An increase of the stimulation 
duration, x e [0, 10], leads to an increase in the desynchro- 
nized transient equally for both, uniform (Figure 10A) and 
optimized, non-uniform stimulation timing (Figure 10B). This 
is indicated by longer intervals of decreased order parameter 
Ri after the stimulation between t on and t 0 {[. For the various 
stimulation lengths, we have f on = 200 and f 0 ff G [205, 217]. 
Beyond x ~ 4.5, the effect of the uniform timing does not 
enhance, while by the optimized protocol it increases further 
until about x = x c ~ 6.6. This value corresponds to the dura- 
tion which assures independence of the stimulation outcome 
on the system's state prior to stimulation (see section 2.4). 
Similarly, small magnitudes of the stimulation intensity 1 < |7| < 
2 yield a prolonged transient for both protocols (Figures 10C.D). 
Beyond |7| = 2, no further increase can be observed for the 
uniform timing (Figure 10C), but the transient after the non- 
uniform protocol continues to grow at least until |7| = 6 
(Figure 10D). 




FIGURE 10 | Influence of the stimulation intensity / and duration x 
on a post-stimulation transient in system (14) with PRC Z = Z ML . 

The four charts visualize the evolution of the order parameter R q for 
x e [0. 10] (A,B), and |/| e [0, 10] (CD). The values of the order 
parameter are encoded in color ranging from 0 (blue) to 1 (red), each 
horizontal strip of a chart corresponds to one time course as shown in 
Figure 9. CR stimulation is applied via m = 4 stimulation sites in the 
stimulation interval [ton, !bff] with t on = 200. The optimized non-uniform 
stimulation protocol is applied in (B,D) and the conventional, uniform 
stimulation timing in (A,C). In (A,B), f 0 ff ranges from f D ff rs 205 to 
f off =a 21 7, and in (C,D), the stimulation interval is constant. In all cases 
the initial phases were drawn randomly from a uniform distribution on 
[0. 2tc]. Other parameters: x = 0.5, N= 240. 



3.3. ROBUSTNESS TO VARIATIONS OF NATURAL FREQUENCIES 

In the above approach we assumed that the neurons are iden- 
tical. In more realistic situations, the parameters of individual 
neurons can vary. In order to test the robustness of the proposed 
resetting technique with respect to parameter changes, which 
break the symmetry of the system, we consider an ensemble (1) 
with non-identical natural frequencies oj ; -,;' = 1, . . . , N, e.g., ran- 
domly chosen from a uniform distribution in [1 — Aw, 1 + Aoo]. 
The results of simulations are shown in Figure 1 1 . We present 
them for the PRC Z = Zml (Figure 1), but qualitatively similar 
results have been obtained for Z = Zh as well. It turns out, that 
a significant prolongation of the post-stimulation desynchroniza- 
tion transient can be observed for a range of Aoo. Indeed, one 
observes a clear difference between the post-stimulation behavior 
of the order parameters R\ and R4 for the suggested CR stim- 
ulation with non-uniform stimulation timing (Figures 11 C,D) 
and that for the conventional CR stimulation with equidistant 
stimulation times (Figures 11A,B). This effect of the optimized 
CR stimulation, however, decreases for a broader distribution 
of the natural frequencies. Nevertheless, our calculations suggest 
that the optimization procedure of CR stimulation can robustly 
improve its desynchronizing impact on neuronal populations 
exhibiting undesired synchronization. 

It is well known that the broadening of the frequency range 
induces a desynchronizing transition (Kuramoto, 1984) such 
that the system with a very broad frequency range does not 




FIGURE 11 | Effect of non-identical frequencies. Time courses of the 
order parameters R q (A,C) and R 4 (B,D) of the neuronal ensemble (14) for a 
range of the frequency detuning Aoo, where the natural frequencies are 
randomly chosen from a uniform distribution in [1 — Atu, 1 + Aoo]. The 
graphical representation is as in Figure 10. The conventional, uniform 
stimulation protocol is applied in (A,B) and the optimized non-uniform 
stimulation timing in (CD). In both cases the initial phases are randomly 
distributed in [0, 2it\. PRC Z = Zml (Figure 1), coupling strength x = 0.5, 
number of oscillators N = 240, number of stimulation sites m= 4, duration 
of single electrode activation x = 10, simulation onset at f on = 200, offset 
at f 0 ff » 217, and stimulation strength / = —10. 
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synchronize even in the absence of stimulation. This transition 
occurs at larger width of the frequency distribution if the coupling 
strength x is increased. In our illustration with coupling strength 
x = 0.5 (Figure 1 1 ) the desynchronization transition and the loss 
of the advantageous effect of the proposed stimulation technique 
take place at approximately the same frequency mismatch. The 
suggested non-uniform CR stimulation is more effective than the 
uniform one for a range of Aw which supports synchronized 
dynamics. 

All results in the sections 3.2 and 3.3 hold within a range of 
the coupling strength, 0 < x < 3.5 for Z = Zml and 0 < x < 5.0 
for Z = Zh (data not shown). In these ranges, with larger x, the 
re-synchronization after the stimulation happens faster, but the 
optimized protocol is still superior to the uniform. This means 
that the length of the transient following the non-uniform stim- 
ulation exceeds the one after the uniform by a similar factor as in 
the presented case for x = 0.5. Moreover, the results are indepen- 
dent of the population size N if it is sufficiently large. In practice 
this is already the case for N > 100, i.e., if each cluster contains 25 
neurons. 

4. DISCUSSION 

A number of pulsatile stimulation techniques have been devel- 
oped which enable to directly shift a synchronized neuronal 
population into a desynchronized state, irrespective of the initial 
state at which the stimulus is delivered (Tass, 2001a,b, 2002a,b). 
However, less favorably, these techniques require careful cal- 
ibration of the stimulation parameters and their continuous 
adaptation to varying model parameters. To overcome this limita- 
tion and provide a desynchronizing stimulation technique which 
is robust and does not require time-consuming or technically 
involved calibration procedures, CR's indirect approach to desyn- 
chronization was developed (Tass, 2003a,b): Inducing a cluster 
state by means of time-shifted phase resetting stimuli delivered to 
different neuronal sub-populations can robustly be achieved and 
does not require relevant calibration (Tass, 2003a,b). The cluster 
states, in turn, are relevant since they lead to long post-stimulus 
desynchronized transients (Lysyansky et al., 2011a), and in the 
presence of STDP (Gerstner et al, 1996; Markram et al., 1997) 
the related decrease of the rate of coincidences induces an anti- 
kindling (Tass and Majtanik, 2006; Hauptmann and Tass, 2007; 
Tass and Popovych, 2012). Neither in preclinical nor in clinical 
studies adverse effects of CR stimulation have been observed (Tass 
et al, 2009, 2012a,b). 

We considered model networks of weakly pulse-coupled neu- 
rons with phase resetting curves and compared them to averaged 
models, where the phase dynamics depends only on the phase dif- 
ferences between the oscillators. Whereas the latter models are 
better analytically tractable and attained a great attention in the 
literature (Strogatz, 2000; Winfree, 2001; Acebron et al., 2005), 
they neglect some important information about stationary states 
of the original systems. In particular, the stationary splay and 
cluster states are not uniform for the pulse-coupled networks, 
contrary to those of the averaged models. These non-uniformly 
distributed cluster states can serve as target states for CR stimula- 
tion in ensembles of pulse-coupled neurons. We have found that 
the optimal stimulation sequence should be slightly non-uniform 



in order to approach the non-uniform cluster state at the end of 
the stimulation. 

We have shown that the phase response curves of the stimu- 
lated neurons determine the phase distribution densities of splay 
and cluster states, which, in turn, influence the timing of the 
stimulation sequences. The proposed non-uniform stimulation 
sequences result in significant improvements of the stimulation 
outcome and lead to several times longer post-stimulation tran- 
sients in comparison to the equally spaced stimulation sequences. 
Intriguingly, modifications of the stimulation timing points of 
only a few percent (e.g., 4%, see Figure 7) actually double or even 
triple the duration of the post-stimulation desynchronization 
transient (Figure 9). The proposed approach takes into account 
and compensates for the interactions among neurons during the 
stimulation. 

We also showed that the discussed stimulation protocol is 
robust with respect to variation of the natural frequencies, stim- 
ulation parameters, and coupling strength. It can lead to a pro- 
longed transient for a range of non-identical frequencies of the 
single oscillators. One can expect that the non-uniform stim- 
ulation technique can be superior to a series of equally timed 
stimulations in more diverse and realistic setups, where trans- 
mission delays and coupling functions are heterogeneous and this 
should be confirmed in further studies. Moreover, since the mech- 
anism of the discussed desynchronizing method is based on the 
phase reset of the neuronal oscillations, which is a universal phe- 
nomenon of the neuronal dynamics (Winfree, 1977; Best, 1979; 
Tass, 1999), and the timing of the optimal stimulation sequence is 
determined by the phase response curves, the presented approach 
is generic and can be applied to other neuronal models and other 
stimulation-induced target states. In particular, to models, which 
employ PRCs of second, or higher order (Oprisan et al., 2004; 
Achuthan and Canavier, 2009). For our purpose, we restricted the 
investigation to the system (1), since it is one of the simplest mod- 
els, which already possesses enough features illustrating our main 
finding, that the optimal stimulation timing is non-uniform. In 
the framework of the considered model one can, however, incor- 
porate the PRCs measured either experimentally or obtained by 
detailed modeling of the globus pallidus and STN (Schultheiss 
et al., 2010; Farries and Wilson, 2012b,a) which are possible tar- 
get regions for CR deep brain stimulation (Popovych and Tass, 
2012; Tass et al., 2012b), and where a change in PRC struc- 
ture might contribute to disease-related changes in synchronous 
activity. Another target for non-invasive acoustic CR neuromod- 
ulation adapted for the treatment of tinnitus (Tass and Popovych, 
2012; Tass et al., 2012a) is the auditory cortex where a phase reset 
can be achieved by different types of auditory stimuli (Brandt, 
1997;Thorne et al., 2011). 

As yet, only slight modifications of the x/m timing of CR 
stimulation (with onset times ti, t\ + x/m, t\ + 2x/m, . . . , t\ + 
(m — l)x/m) have been investigated in the Kuramoto model as 
well as in an ensemble of synaptically coupled FitzHugh-Nagumo 
oscillators modeling spiking neurons in the context of M:N ON- 
OFF CR stimulation, where M cycles with CR are followed by 
N cycles without stimulation (Lysyansky et al., 2011a). For that 
stimulation protocol a x/m timing of CR stimulation was used. 
However, in one variant of this protocol the M-th stimulation 
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cycle was prematurely terminated at the break time f 0 ff, where 
f D ff < x. An optimal choice of t 0 g (e.g., f D ff = 0.44t for a particu- 
lar set of model parameters tested) caused a pronounced increase 
of the desynchronization transient, i.e., an increase of the time 
elapsing to resynchronization by a factor of approx. 2. In contrast, 
inappropriate values of f Q ff induced a decrease (e.g., a halving) of 
the resynchronization time. 

To study post-stimulus desynchronization transients of two 
phase oscillators with time-delayed coupling subject to coinci- 
dent, but phase shifted stimulation the transient time T tr (denned 
as the time it takes a trajectory after stimulus offset to perma- 
nently enter into an g vicinity of the stable phase-locked state) 
was computed (Krachkovskyi et al, 2006). For vanishing time 
delays the phase space of that two-oscillator model is simple 
and the optimal phase shift in the coupling term puts the sys- 
tem's phase difference onto an unstable fixed point (Figure 9a in 
Krachkovskyi et al., 2006). In contrast, for non-vanishing time 
delay the phase space gets considerably more complex, and the 
optimal phase shift puts the system onto a particular point in 
phase space where the system gets trapped by a stable manifold, 
leading to a particularly high transient time T tr . Incorporating 
time delays into the coupling of the model studied in this paper 
will certainly increase the complexity of its phase space. Given 
the results by Ref. (Krachkovskyi et al., 2006), we expect that 
such timed delays may have an impact on the resynchronization 
transient and, hence, on the optimal timing of CR stimulation. 

Another important direction for future analysis comes from 
the fact that biological networks typically comprise neurons of 
different kind. Consequently, our approach has to be extended 
to mixed populations, containing neurons of different type by, 
for example, including inhibitory neurons found in human STN 
(Levesque and Parent, 2005). In this work we considered a sim- 
ple all-to-all coupling topology which provides an easy way to 
obtain a synchronized neuronal dynamics serving as a model 



for the dynamical regimes encountered in Parkinsonian patients 
and monkeys (Nini et al., 1995; Levy et al., 2000). This was 
also motivated by the reported strong functional connectivity 
of tremor-related neurons in STN (Amtage et al., 2009) and 
anatomical intranuclear connectivity as follows from experimen- 
tal and modeling studies (Iwahori, 1978; Kita et al., 1983; Gillies 
and Willshaw, 1998, 2004; Shen and Johnson, 2006), but see 
Ref. (Wilson et al., 2004). The considered weak coupling is sup- 
ported by the observed gradual decay and recovery of pathological 
oscillations at the onset and offset of DBS (Kang and Lowery, 
201 1). CR stimulation has been shown to work for other coupling 
topologies and stimulation setups, e.g., for sensory stimulation 
(Popovych and Tass, 2012; Tass and Popovych, 2012; Tass et al., 
2012a). For further details of the effects of CR stimulation, more 
realistic coupling topologies and sophisticated neuronal models 
as well as connections to other neural populations within the 
basal ganglia and cortical brain areas have to be considered. Also, 
the spatial spread of the stimulation current, 3D effects as well as 
optimized electrode geometries (see e.g., Buhlmann et al., 2011) 
have to be taken into account in future studies. Pursuing such 
studies, our ultimate goal is to come up with CR sequences which 
enable to further minimize the stimulation current for DBS. This 
might contribute to a decrease of the rate of side effects caused 
by stimulation spread to neighboring brain areas. By the same 
token, this might enable considerably smaller geometries of the 
implantable pulse generator due to a significant reduction of 
battery size. 
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APPENDIX Now we estimate F(x) as follows 

A. STATIONARY SOLUTIONS OF THE PHASE DENSITY 

EQUATIONS 
A.1 PULSE COUPLED SYSTEMS 

In this section we show how a stationary solution p (t, cp) = p s (cp) 
for (8) can be found and show that it is unique. It follows from (8) 



F(x)<([ A S ) -x 

\Jo 1 + LZ (-At) x) 



Lx 



that this solution satisfies 2 In (1 + Lux) 

p s (cp) (co + xZ(cp)p s (0)) = C (Al) Hence lim -^ ooF(.x) = -oo. 

Now let us show that (ii) holds. The second derivative of F (x) 

with some constant C. Hence can directl y be calculated and reads 



C 



ft(CP) = a3 + KZ(cp)p s (0)- = 2 ( / gd* ) Fl (x), 



— 3 



The constant C can be determined from the normalization 
condition where 8 <&> *) = (* + l Z (+) *) > °> & MO = ^ Z M0> and 

I M + xZ W p 5 (0) ^ = L (A2) *=(/„ ^) -Jf ^jf 

Combining (Al) and (A2), we obtain the expression (11) for T he sign off" coincides with the sign of Fl We apply the Cauchy- 
p s (cp). Evaluating (11) at cp = 0 and taking into accountZ(O) = 0, Schwarz > s inequa i ity t0 obtain 
the equation for p 5 (0) reads 

, (») - ( r S (f = ^ gi ^ 

\Jo 1 + z; Z MO Pi (°) / 

I 1 II I 3 II /" 2l1 3 2 /" 2lt 

To investigate the resolvability of this equation we consider the — Ir ~ III 2 In ~^\\l 2 J 0 ^ ^ ^ J 0 

function 

Thus, Fi < 0 holds and claim (ii) is proven. 

(r 2n ci\|/ \ We have shown that the equation for the phase den- 

J 0 i -\- X-Z (\|/) x ) shy for the pulse-coupled system (8) has a unique station- 

ary solution. Note that smallness of x was not assumed here, 
which is well-defined and smooth on x e (0, ?) with ? = oo in The first order approximation in x of this density takes the 
case xZ( cp) > 0, and t, = min ~^ ( ^,) otherwise. It is easy to see that form ( 1 2 ) . 

F (0) = J- > 0. Furthermore, we will show that _ „ „..„„-,...„ „„..„. .-„ „„„,-.-..„ ....._.„..„ 

w 2n A.2 SMOOTHLY COUPLED SYSTEMS (3) AND THEIR AVERAGED 



(i) lim x ^ ^ F (x) = — t, < 0 and 



COUNTERPART (4) 



_,, . . ^ „ To unify notations, we write the dynamics for the single neurons 
(ii) t (x) < 0. 

in the limit of large population as 

The conditions (i) and (ii) imply that F (x) has a unique root, 

which corresponds to a unique solution p s (cp) of (8). In order to cp — to + xK (cp, p (t, )) , 

show (i), we treat two cases separately. First, if t, < oo, then the 

smoothness of Z (ty) leads to where 

,. f f 2% d* \ K(m n(t \z{y)fi*G{%)p{t,%)d%, for (3), 

ISlUo l + SZ W x J = °° ^''""iJ^Hfo-Bpfr©*, for(4). 

and, therefore, lim x ^ F (x) = — In the second case, when i; = The phase density equation reads 
oo, we use the following estimates: Since Z(\[r) is smooth and 

Z(0) = 0, one can always find a large enough constant! > 0 such 3tP(f, cp) = — 3<p [p (t, cp) (w + xK(<p, p s (•)))] 

that 

{and its stationary solution p s (cp) fulfills 
\|/, 0 < \l/ < 7T 

2it - \|/, 7T < \|/ < 2ti. C = p s (cp) (oo + xK((p, p s (•))) (A3) 
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with some constant C depending on p s (•). Solving for p s (cp) and 
integration over [0, 2jt] with respect to cp gives an expression 
forC 



C: 



(f 



m + xK($,p s (-)) 



Substitution back into (A3) yields 



Ps(<P) 



Jo 



to + xK(q, p s (■)) 
w + xK{^r,p s (-)) 



2jt 



+ O (k) . 



(A4) 



We see that the stationary solution is generally a x- 
perturbation of the constant state. Without going into 
details, we remark that for small x the right hand side 
of (A4) is a contraction with respect to p s (•) on all sets 
Mp = {feC ([0,2ti]) |/>0, 1/(1! = 1, ||/L<P} of all 
continuous densities on [0, 2tt] which are bounded by some 
fi > j- . This implies that there is a unique stationary solution in 
this set, which can be computed by iteration. 

B. STABILITY OF SYMMETRIC CLUSTER STATES 

In this section, we present a method for calculating characteris- 
tic multipliers of the symmetric cluster states of system (1) (see 
section 2.2.3). To this end, we introduce a discrete time map 
(Bl) which describes the system's evolution in between two sub- 
sequent spikes. The main idea of our analysis is a separation of 
the effects of tangential and transversal perturbations to the clus- 
ter state. This allows us to derive the expressions (B3) and (B5) 
for the characteristic multipliers which can be evaluated easily. 
Finally, we explain how these multipliers behave asymptotically 
as N —y oo. 

Let us consider a time moment fn, when some neuron of the 
ensemble (1) just emitted a spike. We choose its index to be N 
and all other indices in such a way, that the phases are ordered as 
2tt > cpi(f^) > • • • > <PN(tQ ) = 0. Until the next spiking event, 
all neurons advance with the same phase velocity go. Without 
loss of generality, we can assume a) = 1. Then the next neu- 
ron cpi reaches the spiking threshold cpi(fj~) = 2tt at time t\ = 
fn + 2tt — cpi(f^). The other neurons q>k are now located at 
<Pk(h) 

is reset to <p\(tt) 



o 

*PfcC*o~) + 2tt — cpi (fj"). After the spike, the first neuron 



0 and all others to cpj-(tf ) = ^(cp^f^)) = 
\L(<Pk(tg) + 2tt — cpi(^)). Because of the monotonicity of the 
resetting map u,, the order of the phases is preserved as 2tt > 
<P2(*i") > • • • > WN(tf) > <Pi(ff") = 0. We define the firing map 
as F : T N T N , T = [0, 2tt] / {0 ~ 2tt} (the points 0 and 2tc 
are identified as the same point) componentwise by 



Ffc(cp) := [i (<pfc+i + lit - <pi) , 1 < k < N - 1, 



(Bl) 



and .Fjv(<P) := 0. It captures the threshold crossing and spiking of 
the oscillator cpi and shifts the indices k i-> k — 1. This map takes 
the ordered phases after some spiking event 2 it > cpi > - ■ • > 
cp N = 0 and maps them to ordered phases after the successive 



the resetting function |i(cp) to be smooth on T, F is also smooth 
on T N . 

The N-fold iteration of the firing map R := F o ■ ■ ■ o F = F N 
is called the return map. Given an initial state cp with <pi > • • • > 
(Pn, the map -R(cp) returns the state after each neuron has fired 
once. In order to simplify notations, we adopt double indices, 
j := cp(£_i)„ + j, to address the j-th member of the l-th clus- 
ter. A stationary, symmetric m-cluster state is a fixed point ip* = 
F" (ip*) for the n = N/m-th iterate of the firing map, which sat- 
isfies cp^ ■ = tyi with 1 < I < m, 1 < j < n, and ordered cluster 
positions 2it > > • • • > \|/ m = 0. For such a symmetric clus- 
ter state, the time x elapsed between two successive threshold 
crossings of clusters is the same. This time t can be determined 
as the smallest solution of the equation G™~' (0) + x = 2 it. Here 
the function G x (x) := u," (x + x) describes the phase resetting of 
a neuron located at position x + x induced by a cluster spike. 
The value of G'"^ 1 (0) = G T o G x o ■ ■ ■ o G t (0) is the position 
of a cluster initially located at \|r = 0 after m — 1 other clus- 
ters have fired. The positions of other clusters are then given as 

= G™~ £ (0). The linear stability of the cluster state (p* can be 
determined by its characteristic multipliers with respect to R, i.e., 
the eigenvalues of D_R(ip*). In the following we calculate these 
multipliers. Using double indices for F in the same way as for cp, 
we have 

Fe,j (<p*) = W Ok+i + 2ti - th) = to, 

where the index I + 1 is considered modulo m within the range 
from 1 to m. For a perturbed cluster state, the phase of each neu- 
ron can be written as <p£ j = + iu,j- To study the dynamics of 
the perturbations, we introduce the following subspaces 



k=\ 

Vf = p| [r\ kt j = 0, for 1 < j < n] n [r\ e ,„ = 0} , 

kjte 

1 < I < m. Elements of W correspond to perturbations, which 
rule the relative motions of the clusters. Correspondingly, the sub- 
space V[ contain perturbations, which split the l-th cluster. Since 
the perturbations from Vi cannot split any other clusters, we have 



F" (<p* + n) 6 

This implies 

DF" (cp*) Vt C V t 



cp* + W 

cp* + Vi-i ew 



for r) G W, 
for r) G Vi . 



W, and DF" (cp*) W C W. (B2) 



By m applications of the map DF" 
following sequence 



DF" (cp*), we obtain the 



DF" 



DF" 



Hence, for DR (cp*) = (DF" (cp*)) 



DF" 



v«_ m e w = v £ e w. 



we have the skew structure 



spiking event 2ti > -Fi(cp) > • • • > Fn(<P) = 0. Since we assume DR (cp*) V< C >V © Vi and DR (cp*) W C W. This implies that 
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the matrix DR (tp*) has a block triangular form in an appropriate 
basis and the spectrum of eigenvalues of DR ((p*) is the union of 
the spectra, which are restricted to the subspaces W and Vi. The 
eigenvalues belonging to the subspace Vt, can be found by con- 
sidering all possible perturbations in Vi. These perturbations can 
be described as a linear combination of the following basis vectors 



forfc: 
else, 



i < j < i, 



for 1 < i < n — 1. Each element rf corresponds to a different 
way to split the l-th cluster into two subclusters. It is possible 
to show that the multiplier corresponding to rf is given as the 
change of the distance of these two subclusters during one return. 
Calculating the perturbed dynamics step by step (from one cluster 
firing to the next) and linearizing in if leads to 

DR((p*)r\' = \r\' + w 

with some w € W and the eigenvalue X independent on I and i. 



(B3) 



1 = 1 



The perturbations in W can be described by the following basis 
vectors 



1, for k = t, 1 < j < n 
0, else. 



Each element rp corresponds to a perturbation of the £-th cluster, 
which does not split the cluster. Calculations (for brevity, we omit 
the details here) of the dynamics of r| £ show that the eigenvalues 
of DR(ip*) restricted to W are given as the m-th powers of the 
eigenvalues of the matrix 



-8 2 8 2 0 0 

: 0 '■• 0 

-8m : '• 8,„ 
0 0 ••• 0 



(B4) 



This means that tangential multipliers of the symmetric m-cluster 
state are given as X m where X are the solutions of 



det{A - XI) = 0. 



(B5) 



As explained in detail in (Lucken and Yanchuk, 2012), the repeti- 
tive firing |x"(cp) » t}(^, cp) can be approximated by the solution 
cp) of the initial value problem 



(r, cp) = xZ (ft (r, cp)) , with iJ (0, (p) = (p. (B6) 



Furthermore, we have 



dip \m 



.,)) 



m , forZ((p)#0, 
exp (^Z'(cp)) , forZ(<p) = 0. 



(B7) 



This gives an asymptotic method to determine position and sta- 
bility of m-clusters in the limit N — »■ oo by determining the 
quantities 8^ = g| + 2 it — ij/i) for the fixed point ij/ = 

. . . , 4r m ), \Jr i > • • • > \Jr w = 0, of the asymptotic m-cluster 
map F m : T m — > T" ! , which is defined componentwise as 



Fm,l W = * ( — > + 2lt - 1^1 



1 < £ < m. 



(B8) 



This fixed point approximates the stationary, symmetric m- 
cluster <p* in the limit of large N. The stability of cp* can be 
determined by replacing 8^ with 8^ in the expressions (B3) and 
(B4). The asymptotic cluster positions \[r can be determined from 
the smallest solution x of the equation G™ _1 (0) + x = 2jt, where 
G T (<p)=i>(i <p + x). 

C. MODEL EQUATIONS 

The dimensionless Morris-Lecar neuron model (Ermentrout, 
1996; Sato et al., 201 1) is given by the following equations 

V = I - g L (V - V L ) - g K w(V - V K ) 

-gCamoc (V) (V - V Ca ) , 
w = [l\ (V) (Woo — w) , 



with 



1 / , / V - Vi 

moo(V) = - (l + tanh(— — 

1 / /v-V 3 
Woo (V) = - ( 1 + tanh I — — 



1 , / V - V 3 

X (V) = - cosh 

3 \ 2V 4 



Parameters have been chosen according to references 
(Ermentrout, 1996; Sato et al, 2011) as V L = -0.5, V K = -0.7, 
V Ca = 1.0, gL = 0.5, g K = 2, gca = 1.33, V x = -0.01, 
V 2 = 0.15, V 3 =0.1, V 4 = 0.145, 7 = 0.0695 and pi = 0.25. 
We have obtained the PRC Z(cp) = Zml(<p) for perturbations 
V V + AV by direct simulation for AV = 0.0025 and 
setting Zml(<p) = ~Y5V> wnere T IS th e unperturbed period 



of the model and A T is the asymptotic phase lag caused by the 
perturbation. 



Frontiers in Computational Neuroscience 



www.frontiersin.org 



May 2013 | Volume 7 | Article 63 | 17 



